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ABSTRACT 

We study particle dynamics in local two-dimensional simulations of self-gravitating 
accretion discs with a simple cooling law. It is well known that the structure which 
arises in the gaseous component of the disc due to a gravitational instability can 
have a significant effect on the evolution of dust particles. Previous results using 
global simulations indicate that spiral density waves are highly efficient at collecting 
dust particles, creating significant local over-densities which may be able to undergo 
gravitational collapse. We expand on these findings, using a range of cooling times to 
mimic the conditions at a large range of radii within the disc. Here we use the pencil 
code to solve the 2D local shearing sheet equations for gas on a fixed grid together 
with the equations of motion for solids coupled to the gas solely through aerodynamic 
drag force. We find that spiral density waves can create significant enhancements in 
the surface de nsity of solids , equivalent to 1 — 10cm sized particles in a disc following 
the profiles of Clarke (2009) around a ~ 1M star, causing it to reach concentrations 
several orders of magnitude larger than the particles mean surface density. We also 
study the velocity dispersion of the particles, finding that the spiral structure can result 
in the particle velocities becoming highly ordered, having a narrow velocity dispersion. 
This implies low relative velocities between particles, which in turn suggests that 
collisions are typically low energy, lessening the likelihood of grain destruction. Both 
these findings suggest that the density waves that arise due to gravitational instabilities 
in the early stages of star formation provide excellent sites for the formation of large, 
planetesimal-sized objects. 
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1 INTRODUCTION 

There are currently two models for the formation of gas gi- 
ant planets in circumstellar discs. The most wi dely accepted 
is the core accretion model (Pollack et al. 199^). Here, a core 
of solid material grows via a series of collisions until it be- 
comes m assive enough to ac crete a gaseous envelope from 
the disc l|Pollack et al.l lT996). This must occur before the 
disc is depleted of gas, a proc ess which is obser vationally 
estimated to take ~ 10 7 years (|Haisch et al.ll200lh . 

A major area of uncertainty in the core accretion 
model lies in the growth of objects from small dust grains 
to kilometre-sized planetesimals, which form the building 
blocks of planet cores. Many numerical models of the core 
accretion process assume that a large population of plan- 
etesimals has already formed. 

The dynamics of the smaller particles that ultimately 
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grow to form these kilometre-sized planetesimals is, however, 
governed by the drag force that arises from the velocity dif- 
ference between the solid particles and the surrounding gas. 

Within a few scale heights of the disc mid-plane, the 
pressure gradient within the disc tends to be negative. This 
causes the gas to orbit with sub-Keplerian velocities. The 
dust is not affected by the gas pressure gradient and orbits 
at Keplerian velocities. The drag force on the dust particles 
that arises from this velocity difference results in the solids 
losing angular momentum to the disc and d rifting inward at 
a rate that depends on the particles' size jWeidenschillind 
ll977tlNakagawa et al.lfl986l ). For very small grain sizes, the 
dust is tightly coupled to the gas in the disc and the ra- 
dial drift velocities are small. For very large objects, the 
solids are decoupled from the gas, move in approximately 
Keplerian orbits and again have very small drift velocities. 
Particles in the intermediate size range can, however, have 
large drift velocities. Although the exact size range depends 
on the local properties of the disc, drift velocities can ex- 
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ceed 10 3 cm/s for obje cts with sizes between 1cm and lm 
<|Weidenschillingill977l ). In order to prevent dust particles 
from spiralling inward into the central star before it decou- 
ples from the gas, solid material must rapidly grow from 
small centimetre- s ized g rains to large decametre sized ob- 
jects. lLaibe et al.l (|2012T ) do however suggest that there may 
be surface density and temperature profiles for which parti- 
cles may survive this inward migration. 

The alternative disc instability model requires that 
planets be able to form via direct gravitational collapse of 
the gaseous part of t he disc, an d does not require the pres- 
ence of a solid core (|Bosslll998l ). Gravitational instabilities 
set in, for a razor-thin disc, when the sound speed c s , Ke- 
plerian rota tion frequency , SI and the surface density of the 
disc satisfy <|Toomrelll964h 



C 3 Si 

ttGE 



< 1. 



(1) 



This requires that the disc be relatively massive compared 
to its parent star. In the event that a disc is susceptible 
to such instabilities, one of two outcomes may occur. If 
the cooling time is long, the disc will settle into a quasi- 
steady state, where the coo ling balances t he heating gener- 
ated by gravitoturbulence ijGammi el l200ll ). For short cool- 
ing times the disc may fragment, forming gas giant plan- 
ets. The critical cooling time below which frag mentation 
occurs is commonly t aken to be t c , cr it = 3S1 -1 ([Gammiel 
l200ll ; [Rice et al]|2003l ) however recent studies suggest that 
this threshold may not be fully converged, with recent 
high resolution simulations suugesting that the critical coo l- 
ing time, t c>cr i t may exceed 10fl -1 (|Meru fc Bate! [20 111 ). 
It has, however, been suggested that this result i s nu- 
merical JPaardekooper et al. 2011 1_ Lodato fe Clarke] l201ll ; 



iRice et al.ll2012l ). |Paardekooperl (|2012T T do, however, suggest 
that there may be an intermediate range of cooling times 
for which fragmentation may be stochastic, observing frag- 
mentation in some simulations with cooling times as high as 
to = 20fi _1 . 

Although very few Class II protostars are observed to 
have sufficiently m assive discs to be suscepti ble to gravita- 
tional instabilities (|Beckwith fc Sargend[l99"lh . observations 
indicate that during the Class and Class I phases, mas- 
sive discs may be much more common a r ound protostars 
jRodriguez et all 120051 ; lEisner et all 120051 ; [g reaves fc Rice! 
120101 ). However, even though a massive circumstellar disc is 
thought to be present at early times in the star formation 
process l|Machida fc Matsumotdl201ll ). it is not clear that 
such discs can cool quickly enough to directly form giant 
planets via fragmentation. 

Even if the cooling times present within discs are too 
long to allow giant planets to form via fragmentation, it 
is expected that circumstellar discs are self-gravitating i n 
their early stages |Lin fc Pringle.|["l987l ; lLin fc Prin glc 1990). 
If this is the case, these instabilities take the form of non- 
axisymmetric spiral structures, which have been shown to be 
highly efficient at transporting angular m omentum outward 
jLodato fc Ricdl2004l ; iForgan et al]|201ll ) . 

It has been proposed that these spiral waves may be 
highly effective at trapping the solids in the disc. The gas 
pressure gradient changes from positive to negative across 
these spiral wave structures, which results in sub-Keplerian 
velocities on one side of the wave, and super-Keplerian on 



the other. The drag force then causes dust grains to drift 
toward the density/pressure maxima at the centre of the 
waves. This can eliminate two problems currently challeng- 
ing theories of planet formation. Firstly, by trapping large 
amounts of solid material within a spiral arm the local den- 
sity of solids is much higher than the average within the 
disc, leading to the faster creation of planetesimals. Trap- 
ping particles in these local pressure maxima can also shield 
growing objects from the rapid inward drift described previ- 
ously that can potentially stop la rge objects form ing before 
they drift into the central object. IRice et alJ (|2004h showed, 
using global simulations, that the surface density of certain 
particle sizes can be e nhanced by a fact or of over 100 in the 
spiral wave structure. IRice et alJ (|2006l ) estimate that such 
an augmentation in the surface density of solids will lead 
to the creation of k m-scale planet esimals through the self- 
gravity of the solids. Idarke j (120091) howeve r not ed that the 
simulations carried out in IRice et al l (120041 1 and[l ice et al.l 
(2006) assumed a constant cooling time throughout the disc, 
and that the conditions required to achieve such high solid 
concentrations may be restricted to very large orbital radii 
in discs that incorporate a more realistic cooling. 

Similar concentrations of particles are also seen in 
the presence of other density enhancements in the gas, 
such as those cause d by magneto-rotati onal turbulence 
Jjohansen et alj|2006l), as well as vortices (|Godon fc Livid 
|2000| ; iKlahr fc Bodenheimeij 120031 ). gaps in the gas that 
form due to the presence of a massive planet are also ob- 
served to create density enhance ments in the solid compo- 
nent of the disc at their edges JPaardekooper fc Mellemal 
120061 ; iFouchet et aHl2007l ; lAvliffe et alj|2012l ). Increased par- 
ticle concentrations are thought to aid planetesimal forma- 
tion in two ways , the enhanced collision rate displayed by 
I Rice et all (|2004l ) will aid planetesimal formation (provided 
collisions are constructive). Secondly, the presence of high 
dust densities may actually lead to the formation of large 
planet esimals through direct gravitational collaps e of the 
solids l|Rice et al. 2006; [Johansen et al. I l2007l , l201ll ). 

The ultimate goal of the present paper is an improved 
understanding of grain growth and planetesimal formation 
at early evolutionary stages when the influence of the disc's 
self-gra vity is non-negligable. Spe cifically, motivated by the 
work of lClarke and Lodatol l|2009i ). we study dynamical be- 
haviour of particles embedded in a self-gravitating disc using 
a local shearing-sheet app roximation. As mentioned above, 
IClarke and Lodatd (2009) concluded, based on analytical es- 
timations, that in realistic self-gravitating discs, particle ag- 
glomeration will be restricted to the outer regions of the disc 
(at radii of several tens of AU). At such large separations, 
the cooling time-scale is relatively short and hence density 
enhancements in spiral features are strong enough to con- 
centrate particles on a local dynamical (orbital) time scale - 
the lifetime of such features. Here, we go beyond these esti- 
mates and study in greater detail, via numerical simulations, 
the trapping capability of spiral waves for a wide range of 
disc cooling times and particle stopping times. This allows 
us to see how the solids will respond to the presence of spiral 
density waves at various radii within the disc and also un- 
derstand what minimum values of gas density enhancements 
are required for appreciable particle concentration to occur 
is spiral features. 

This paper is arranged as follows. In Section 2 we 
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outline the model we use in our simulations. In Section 3 we 
discuss the evolution of the gas and dust particles. Finally, 
we draw our conclusions in Section 4. 



2 MODEL 

To investigate the dynamics of solid particles embedded in 
self-gravitating proto-planetary discs, we solve the 2D lo- 
cal shearing sheet equation s for gas o n a fi xed grid, includ- 
ing disc self-gravity as in iGammi together with 
the equations of motion of solids coupled to the gas solely 
through aerodynamic drag force. As a main numerical tool, 
we employ the pencil codeQ. The pencil code is a sixth 
orde r spatial and third o rder temporal finite difference code 
fsee lBrandenburd (|2003h for full details). The pencil code 
treat s solids as numerical super-particles (|johansen et al.l 
120061 ). 

To use the local shearing sheet model, we must make 
two simplifying assumptions, first that the disc is cool, and 
therefore thin ( H/r ~ c 3 / ( Q,r) <C 0.1), and that the disc is 
razor-thin, as in iGammiel (|200ll ). In the shearing sheet ap- 
proximation, disc dynamics is modelled in the local Carte- 
sian coordinate frame centred at some arbitrary radius, ro, 
from the central object and rotating with the disc's angu- 
lar frequency, Q, at this radius. In this local frame, the x- 
axis points radially away from the central object, the y-axis 
points in the azimuthal direction of the disc's differential ro- 
tation, which in turn manifests itself as an azimuthal parallel 
flow characterised by a linear shear, q, of background veloc- 
ity, uo = (0, — qQ,x). Our simulation domain spans the re- 
gion —L x /2 < x ^ L x /2, —Ly/2 ^ y < L y /2. We adopt the 
stand ard shearing-sheet boundary conditions (jHawlev et al.l 
1 19951 ). namely for any variable / we have 

f(x,y,t) = f(x,y + L,t) (2) 

f(x,y,t) = f(x + L,y-qSlLt,t). (3) 

These boundary conditions apply to all grid variables except 
the azimuthal component of the velocity, u y , which must be 
adjusted to account for the relative shear between neigh- 
bouring boxes. The azimuthal velocity on the radial bound- 
ary is, 

v y (x,y,t) = v y (x + L,y — qflLt, t) + qflL. (4) 

The shear parameter q = 1.5 for Keplerian rotation profile 
adopted in this paper. 

2.1 Gas Density 

In this local model, the continuity equation for the gas sur- 
face density E is 

^ + V.(Eu)-gnx^-/ D (E) = (5) 

where u(u x , u y ) is the gas velocity relative to the back- 
ground Keplerian flow, which in our domain manifests it- 
self as an azimuthal flow with linear shear of velocity uo = 

1 Sec http://code.google.eom/p/pencil-code/ 



(0, — gOx). Due to the high-order numerical scheme of the 
pencil code we also include a diffusion term, fo , to ensure 
numerical stability and capture the effects of shocks, 

/D=<i?(V 2 E + Vln<i5-VE). (6) 

Here the quantity £d is defined as, 

( D = D sh (max[(-V.u) + ])Ax 2 (7) 

where D s h is a constant defining th e strength of shoc k dif- 
fusion as outlined in Appendix B of lLvra et alj (2008). 

2.2 Gas Velocity 

The equation of motion for the gas velocity u relative to the 
background Keplerian flow uo takes the form 

c* u / „ x „ du VP „„„ _ 

— h (u ■ V)u — qilx— — — 2S2z x u + q\lu x y 

at ay E 

+ 2fiAux- VV + f^u), (8) 

where P is the two-dimensional pressure and xj) is the gravi- 
tational potential of the sheet due to the gas surface density 
perturbation E — Eo. The left hand side of equation [8] in- 
cludes terms from the velocity field u and the Keplerian 
flow. The first term on the right hand side of the equation 
is the force due to the pressure gradient. The second and 
third terms represent the Coriolis effect/shear induced by 
the choice of coordinate system. The fourth term mimics 
a global radial pressure gradient, which reduces the orbital 
speed of the gas by Av and is responsible for the inward 
migration/drift of solids in an unperturbed disc. The main 
purpose for the inclusion of this term here is to see how the 
radial drift affects the concentration of particles by spiral 
density waves. The fifth term represents the effect of the 
gravitational potential of the disc. Finally we include an ex- 
plicit viscosity term 

f„ =!/(V 2 u + iw ■ u + 2S ■ Vln E) 

+ C4V(V.u) + (Vln E + Vln C,)V.u], (9) 

which contains both Navier-Stokes viscosity and a bulk vis- 
cosity for resolving shocks. Here S is the traceless rate-of- 
strain tensor 

I ( dm duj 2 \ 
^ = 2{dx- + dx-- 3 5 - V - U J (10) 

and is the shock viscosity coefficient analogous to the 
shock diffusion coefficient defined in equation [7] with v s h = 

D 3h 

2.3 Self-gravity 

The gravitational potential of the gas is found by inverting 
the Poisson equation for a 2D (razor thin) disc 

V 2 t/> = 47rG(E-Eo)5(2) (11) 

using the Fast Fourier Tr ansform method outli ned in the 
supplementary material of lJohansen et all (|2007l ). Here, the 
surface density is Fourier transformed from the (x, y) plane 
to the (k x ,k y ) plane withou t the intermedi ate co-ordinate 
transformation performed bv lGammiel (|200ll ). In this case a 
standard FFT is adapted to allow for the fact that the radial 
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wavenumber k x of each spatial fourier harmonic depends on 
time as k x (t) — k x (0) + qQk y t in orde r to satisfy the shearing 
sheet boundary conditions (see also iMamatsashvili fc Ricei 
2009). Note that since the gravitational potential is associ- 
ated with density fluctuations, only the perturbed surface 
density S — So enters equation 1111 



2.4 Entropy 

The pencil code uses entropy, s, as its main thermody- 
namic variable , rather than internal energy, U, as used by 
iGammid |200ll i. The equation for entropy evolution is 



Ec 2 



|-^| + (u.V) S =i(2E^-^-^ 



(12) 

where the first term on the right hand side is the viscous 
heating term and the second term is an explicit cooling term. 
Here we assume the cooling time t c to be constant through- 
out the sheet and take its value to be sufficiently large that 
the disc does not fragment and achieves a quasi-steady state. 
The final term on the right hand side is a shock dissipation 
term analogous to that outlined for density. 



2.5 Dust Particles 

The dust particles are treated as a number of massless nu- 
merical tes t-particles, usin g the implementation for dust 
particles of IJohansen et al] (2007). In this implementation 
dust particles have positions x = (x p ,y p ) on the grid and 
velocities v = (y x , v y ) relative to the unperturbed Keplerian 
rotation velocity uo = (0, — qQx p ) in the local frame. These 
are evolved as 



dx 
di 



— = v - qQ,x p y 



dv 1 

— = — 20z x v + qQ.v x y H (u - 

dt r f 



(13) 



(14) 



wher e r/ is the friction time of the particle jjohansen et al.l 
2006). The first two terms in equation Q3] represent the Cori- 
olis Force and the non-inertial force due to shear. The final 
term describes the drag force exerted by the gas on the par- 
ticles which arises from the velocity difference between the 
two. Unlike the gas the particles do not feel the pressure 
force. The motivation for this paper is to investigate where, 
in a typical self-gravitating disc, the solid particles are most 
strongly influenced by the self-gravitating structures in the 
gas disc. Additionally, we wish to investigate what range 
of particle sizes are most strongly affected, without mak- 
ing any assumption as to how particle sizes and masses are 
distributed. Consequently, we use massless 'test' particles. 
Therefore, there is no term modelling the back-reaction of 
the drag force in equation[8] nor is the self-gravity of the par- 
ticles considered in equation 1141 We do acknowledge, how- 
ever, that if the local solid density became comparable to 
that of the gas, these terms would become important. 

The drag force on the particles from the gas is calculated 
by interpolating the gas velocity field to the position of the 
particle, using the second order spline interpola tion outlined 
in Appendix A of lYoudin and Johansenl l|2007l ). 



2.6 Units and Initial Conditions 

We normalise our parameters by setting c s ,o = ^ = Eo = 1. 
The time and velocity units are [t] = and [v] = c s ,o, re- 
sulting in the orbital period, T — 2n. The unit of length 
is the scale height, [I] = H — c s ,o/Sl. We initialise the 
Toomre-Q parameter to be 1 throughout the sheet. This 
sets the gravitational constant G = 7r -1 . The surface den- 
sity of gas is initially set to be unity everywhere in the 
sheet. The box is of length L = 80GE/O 2 and is divided 
into a grid of 1024 2 cells. This choice of units sets the 
sheet length L = 80H/tt ~ 25H. This is a much larger 
value for L than that typically take n for simulations inves- 
tigat ing MRI driven turbulence (e.g. IJohansen et al] (2006, 
l201lf )). however is t ypical for those investigating instabilities 
due to self-gravity. ijGam mic 2001) estimates that structure 
in self-gravitating discs typically arises on length scales of 
~ L/5 = 16GE/Q 2 . Typical proto-planetary dis cs are es- 
timated to have aspe ct ratios H/r = 0.05 — 0.1 (|Gammiel 
l200ll; lArmitageM201ll). T his is supported by observations 
by I Andrews et all (|20l"(ih . who find that typical discs have 
scale heights ranging from 2-20AU at radii of 100AU, giving 
aspect ratios ranging from 0.02 to 0.2. A shearing sheet cen- 
tred at ro = 20AU with a scale height of ~ 1AU will there- 
fore span the radial range ~ 10 — 30AU. It is worth noting 
that the cooling time, t c , which we have assumed to be con- 
stant througho ut the s heet, in reality is t c = i c (E, U, Q) as 
described by IGammid l|200ll ). However the use of constant 
cooling time over a sheet of this size allows us to infer the 
general behaviour of the dust particles at a given location 
within the disc. 

The gas velocity field is initially perturbed by 
some small random fluctuations with the rms amplitude 
V < <5v 2 > = 10 -3 . We take the viscosity and diffusion 
coefficients to be v — 10 -2 and v S h = D s h = 5.0. We 
use 5 x 10 particles, split evenly b etween five f riction 
tim es, T f = [0 01,0.1 ,1,10, lOOin- 1 . iBai fc Stone] (|2010D 
and lLaibe et al.1 (I2012T ) have shown there is a spatial reso- 
lution criteria which applies to coupled dust and gas sim- 
ulations such as those outlined above. For the dust par- 
ticles to be properly resolved, the grid spacing must sat- 
isfy Ax < c 3 Tf. For the chosen set of parameters we have 
Ax/c s ~ 0.07, they satisfying this condition for all but the 
Tf = 0.01S7 -1 particles. The required increase in resolution 
to satisfy this for the 77 = 0.01J7 -1 particles is beyond our 
computational limits. The effects of these particles poten- 
tially being under-resolved are discussed in section 3. Each 
particle species with a fixed radius is initially given zero ve- 
locity relative to the Keplerian flow, v(t = 0) = (0,0), and 
is distributed spatially uniformly throughout the sheet. 



3 RESULTS 

3.1 Gas Evolution 

The evolution of the gaseous component of the disc is 
in good agreement with that observed in analogous stud- 
ies based on the shearing sheet formalism ( Gammic 1200 lj ; 



iJohnson fc GammielbOOa - ] iMamatsashvili fc Ricel200g| ). The 
small initial velocity fluctuations grow and develop into non- 
linear fluctuations in velocity, surface density and potential. 
Shocks then develop which proceed to heat the gas, while 
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Figure 1. Logarithmic surface density of the gas after 40 orbits 
in the t c = lOQ^ 1 run. At this stage the disc is already in a 
quasi-steady state 



the cooling works to reduce the entropy of the gas. Den- 
sity structures develop which are sheared out by differential 
rotation. These density structures tend to take on a trail- 
ing rupture, which leads to a finite shear stress parameter 
a l|Gammiel l200ll ). After a few orbits, the heating due to 
shocks is balanced by the cooling term and the thermal en- 
ergy of the sheet settles into a quasi-steady self-regulated 
state, where the thermal, kinetic and gravitational energies 
of the disc are approximately constant with time. In this 
state, on the surface density field we clearly see elongated 
trailing surface density features, or density waves (Figure 
[1} , whose particle trapping capability is our main subject of 
study. 

The magnitude of these nonlinear density enhance- 
ments is determined by the cooling time - the smaller the 
cooling time, the l a rger the amplitu de of density waves 
i|Cossins et al.l 120091 ; I Rice et al. Il201ll ). To analyse the ef- 
ficiency of particle trapping by density wave structures at a 
range of radii within the disc, we consider several values of 
t c — 10, 20, 40, 80 and 160f2 -1 , bearing in mind that in discs 
realistic (radiative) effective cooling time - scales tend to de- 
crease with increasing radius l|Clarkeil2009l ; lRice fc Armitaed 
l2009h . 



3.2 Particle Concentration 

The particles, which initially are given random positions 
and zero velocities (relative to the Keplerian flow), are not 
evolved until the gas has undergone its initial burst phase 
and settled into a quasi-steady state. Once the gas has 
reached this phase, we release the dust particles. For the 
t c = 10n _1 and 20tt~ 1 runs, the particles are introduced at 
a time t par = WT, while for the t c = 40f2 _1 , 80fi _1 and 
160r2 _1 runs at t par = 20T, t par = 30T, and t par = 50T re- 
spectively. In each case we evolve the system for a further 30 
orbits, until the particles have also reached a quasi-steady 
state and come into a dynamical equilibrium with the gas 
- the majority of particles are trapped in density wave fea- 



tures, but particles aggregated within each density enhance- 
ment disperse as the the latter gets sheared out by the disc's 
differential rotation. These particles are then captured in a 
new density structure. 

Upon release, the particles are drawn to local pressure 
maxima associated with the density waves in the gas. Fig- 
ure [2] shows the logarithmic surface densities of the dust 
grains with different friction times in a quasi-steady state 
at t — 40T in the t c = 10£1 — 1 run. The corresponding gas 
surface density at the same time is plotted in Figure [1] Com- 
paring Figures [T] and [5] we see a clear correlation between 
the density enhancements in the gas and over-densities in 
the particles. 

The degree to which the presence of spiral density waves 
affect the particle dynamics and the magnitude of their con- 
centration depends on the friction time of the particles, as 
evident from Figures 2 and 3. The first three panels in Figure 
2 show the surface densities of the smaller particles with cor- 
respondingly smaller friction times (77 = [0.01, 0.1, 1.0]f2 -1 
respectively). In this case, there is a high degree of correla- 
tion between the structures in the gas surface density and 
the arrangement of the particles, with the smallest parti- 
cles (with 77 ^ O.lfi -1 , top two panels) tightly mapping the 
overall gas structure but the local enhancements of parti- 
cle density are not exceptionally large, with typical density 
enhancements of no more than 10 times the mean particle 
density. The intermediate-sized particles (with 77 = l.QQ~ , 
middle left panel) are primarily concentrated at the crests of 
spiral waves, tracing out the locations of the pressure max- 
ima and reaching there the largest values of density, up to a 
factor of ~ 10 3 times the mean particle density. By contrast, 
there is much less correlation between the surface densities 
of the gas and the larger particles (with 77 l.Ofi -1 , mid- 
dle right and bottom left panels) and therefore the particle 
over-densities are also lower. In the latter case, the mean 
Keplerian motion of the particles tends to dominate over 
the action of the drag force, reducing the effect of the lat- 
ter to a small perturbation of the background motion. The 
specific behaviour of particles with various friction times 
(sizes) in the density field of spiral waves found here agrees 
well w ith that in global disc simulations of particle dynam- 
ics by IRice et al.l 1 20041 ). In this connection, we would like 
to mention that a similar high degree of concentration of 
intermediate-sized (~ lm and 77 ~ 1.0O -1 ) dust particles 
also occurs in gas over-densities produced by MRI-driven 
turbulence in discs l|johansen et al.1120061 . 120071 ). 

Thus accumulation of the particles within spiral waves 
leads to an enhancement in the local surface density of the 
solids within the disc. Figure [3] shows how the maximum 
value of particle density, relative to the mean particle den- 
sity, inside the domain varies with time for a range of particle 
sizes (friction times) in the t c = lOfi -1 simulation. Particles 
with very short (17 ^ 0. 1£1 — 1 ) and very long (77 ^ 10S7 -1 ) 
friction times exhibit small density enhancements, with the 
maximum density of particles typically no more than a few 
times the mean (initial) density. However, particles that oc- 
cupy an intermediate size range with 77 = 0.1 — 10f2 -1 ex- 
hibit orders of magnitude larger growth in gas over-densities. 
For example, particles with 77 = 0.1S7 -1 and 77 = lOQ.^ 1 
experience an increase in concentration approximately 10 
times larger than that of particles in the above mentioned 
first two size ranges and around 20-30 times that of the 
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Figure 2. Logarithmic surface density of the particles for the t c = 10f!~ 1 run at t = 40T, when the disc is already in a quasi-steady 
state, 30 orbital periods after the drag force between the gas and particles has been turned on. The surface density of gas at this time 
is shown in Figure 1. The first five panels show the surface density of particles with friction times of tj = [0.01, 0.1, 1.0, 10.0, 100. 0]fi -1 , 
respectively. The final panel shows the total surface density of the dust summed over all particle species. We see that intermediate-sized 
particles with Tf ~ 1.0r2 _1 arc captured most effectively in density structures. 



Planetesimal Formation In Self- Gravitating Discs 7 



initial particle density throughout the sheet. As mentioned 
above, particles characterised by 77 = 1.0f2 _1 are remark- 
able as they exhibit the highest degree of concentration in 
density wave crests/maxima, with the local density of the 
particles typically reaching levels a factor of ~ 10 2 times 
that of the mean value, or initial particle density. This range 
of particle density enhancements is in good agreement with 
analytic theory, which predicts that the effects of the aero- 
dynamic drag force are most significant for 77 = 1.0fl~ 
particles. This consistency with analytic theory is reassuring 
when considering the potential non-convergence issue with 
the smallest (77 = O.OIQ." 1 ) particles, as combined with the 
lack of obvious numerical problems in Figure 2, suggest that 
our results are correct to a good degree of accuracy, even for 
the insufficiently resolved particles. 

In this work we have used massless 'test' particles so as 
to focus mainly on where the self-gravitating structures in 
the gas disc can have the most significant affect on the solid 
particles. In a typical star forming cloud, the solid particles 
make up ~ 1 % of the mass. This can, however, be dis- 
tributed over a wide range of particle sizes and so it is not 
clear how much mass will be in particles that are strongly 
influenced by the gas overdensities (i.e., in particles with fric- 
tion times 0.ir2 _1 < 77 < 10fi _1 ). Given that the density of 
these particles can be enhanced by a factor of 100 — 1000, if 
~ 10 % of the solid particle mass is in this size range, there 
will be regions in the disc where the gas and solid densi- 
ties are comparable and we would then need to include the 
backreaction of the solid particles on the gas and also the 
self-gravity of the solid particles to determine their subse- 
quent evolution. 

The enhancement in the local density of solids is also 
dependent on the gas cooling time, t c . As outlined in Section 
13.11 the amplitude of spiral wave structure is determined by 
the value of cooling time, with smaller cooling times result- 
ing in larger pressure gradients and higher surface densities 
in the spiral waves. The higher pressure gradients due to 
short cooling times tend to be more efficient at trapping the 
particles. Figure [4] shows the maximum of net surface den- 
sity of all particles over time for a range of cooling times. 
Although in this plot we show the net particle density, at 
any given value of cooling time the maximum value of den- 
sity is still due to the particles with 77 = 1.0ft , similar 
to the above situation with t c = 10fi _1 . Here the amount 
of particle accumulation that occurs appears to have a de- 
pendance on t c , with the longest values studied (t c = 80f2 _1 
and t c — lGOfl^ 1 ) showing less concentration than that for 
the shorter cooling times. This decrease is due to the shal- 
lower gradients in both pressure and potential associated 
with the longer cooling times. Although longer cooling times 
are not considered here, we expect this trend to continue, 
with longer cooling times exhibiting less concentration. The 
particles also reach a quasi-steady state much quicker in the 
simulations with a shorter cooling time. As the cooling time 
gets longer, the particles tend to form long filaments of ap- 
proximately uniform density along the spiral waves in the 
gas, rather than tightly crowding in density wave crests. 

It is instructive at this point to relate the parameters 
used in our simulations to realistic models of self-gravitating 
proto-planetary discs. By taking the analytically derived 
profile for the total stresses a from the appendix of IClarkel 
l|2009l ) and finding with it an effective cooling time using 
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Figure 3. Maximum particle density, relative to the mean, as a 
function of time for a range of friction times in the t c = 10f2 _1 
simulation. Particles with 77 = 1.0f2 — 1 exhibit the highest degree 
of concentration in gas over-densities. 
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Figure 4. Maximum particle density, relative to the mean, as a 
function of time for a range of cooling times. Particles are some- 
what more effectively concentrated in spiral density waves at short 
cooling times, with the t c = 80f2~ 1 and the t c = 160f2 -1 sim- 
ulations showing a decrease in particle concentration than the 
shorter times considered. 



the expression given in iGammid |200ll ). we can estimate 
the range of radii in the disc corresponding to the effec- 
tive cooling times considered here. Figure [5] shows how such 
obtained t c varies with radius within the disc, for given stel- 
lar accretion rates M. From Fig. [5] we see that the range 
of cooling times considered here spans a radial interval 
20 — 60AU in the disc. This is actually the range of radii 
at which disc self-gravity is appreciable (Q ~ 1) and main- 
tains the disc in a quasi-steady state. At larger distances 
70AU) cooling is too fast, causing the disc to fragment, 
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Figure 5. Analytically derived profile for t c as a function of ra- 
dius for different mass accretion rates. The effective cooling times 
we consider here correspond to the radial range ~ 20 — 60AU. 



whereas at smaller distances 10 AU) cooling is very inef- 
ficient and so the self-gravitating density perturbations will 
be negligible and this region is lik ely to be dominated by 
the magneto-rotational instabilit y |Balbus fc Hawle"vlll99ll ; 
IZhu. Hartmann fc Gammie!l2009r i. Once a radius within the 
disc is specified, we can also determine the actual values 
of physical variables (gas density, sound speed, etc.) in our 
simulations and relate the dimensionless friction times of 
particles used here to their real physical sizes. 

The force acting on a dust particle of radius a traveling 
at a velocit y u relative to the surrounding gas o f density p 
is given by |Whipple!l 19721 IWeidenschillindll977r i, 

1 



-Gn^a pu u, 



(15) 



where u = |u|, u — u/u and the drag coefficient Cd is de- 
fined as, 



8 c, 
3 u 



a < 



9A 



C D = { 



24R7 1 R e < 1 



(16) 



24R: 



0.44 



1 < i? e < 800 



R e > 800 



for low Mach numbers. The drag regime with a < 9A/4 is 
generally called the Epstein regime, whereas the other three 
regimes define the Stokes drag. In expressions 1161 R e is the 
Reynolds number characterizing gas flow in the vicinity of 
a dust particle, which will be defined below, and A is the 
mean free path of gas molecules. Assuming the gas to be 
made mainly of molecular hydrogen, for the mean free path 
we get 



A = 



rriH., 



4 x 10" 9 



pA p 
The Reynolds number is given by, 

Ipau 



Re 



'I 



(17) 



(18) 



where r\ = pv is the gas viscosity. For collisional viscosity, 
we have 

v = *£- d9) 

The Reynolds number can therefore be expressed as 

* -«($)(£). w 

Once the drag force on the particles is known, the corre- 
sponding friction time of the particles (as used in equation 
I14[) can be calculated, 



(21) 



where m p = |7ra 3 p p is the mass of a particle of internal 
density p p . 

By following this prescription, we estimate t hat for 
a disc with the surface density profiles derived in IClarkel 
(2009), at the innermost radius considered (~ 20AU), the 
range of particle sizes considered spans Imm-lOm, with the 
tj — 1 particles corresponding to objects with a physical 
size of 10cm. At an outer radii of ~ 60AU, the size range 
considered is halved, spanning 0.3mm-3m, with 77 = 1 par- 
ticles corresponding to 3cm size objects. We note here that 
the inner radii tend to be more effective at collecting 'larger' 
particles, this allows for the possibility of small objects form- 
ing at large radii and growing as they drift inward due to 
aerodynamic drag. 

3.3 Particle Velocities 

Although large concentrations of dust and ice is vital for the 
growth of planetesimals, the velocities of the particles must 
also be considered in order to determine the evolution of the 
solid component of the disc. If the relative velocities of the 
particles are too large, collisions become too energetic, and 
tend to result in large objects being broken apart, rather 
than growing. Figure [6] shows the distribution of velocities 
for each particle size at the end of the t c — lOfi -1 simula- 
tion. The peak of the distribution lies at higher velocities for 
smaller grain sizes, whilst the width of the distribution tends 
to also increase. In the case of the 77 = l.Ofi -1 particles, the 
width of the distribution is unusually narrow. This appears 
to be a result of the particles being tightly concentrated at 
the centre of density waves, causing their velocities to be- 
come highly ordered, resulting in the majority of particles 
having similar velocities. This keeps the relative velocities 
between particles low, improving the likelihood of collisional 
growth occurring. 

By contrast, particles with small 77 ^ 1.0£7 — 1 are 
strongly coupled to the gas and practically repeat the over- 
all turbulent motion of gas outside spiral arms, hence the 
larger dispersion in the particle velocity. On the other hand, 
particles with large 77 ^ 10. Oil -1 are very loosely coupled 
to the gas and therefore most of them have small velocity 
dispersion, i.e., they are not effectively accelerated by gas. In 
Fig. 6, this corresponds to the large peak in particle number 
near zero velocity for large friction times. 

As we move radially inward within the disc, we see the 
velocity dispersion for each particle size narrows as the cool- 
ing time increases. Figure [7] shows the velocity dispersion 
of the 77 = l.Ofi -1 particles for a range of cooling times. 



Planetesimal Formation In Self- Gravitating Discs 9 



, x 10 



x = 0.01Q 



= o.m _ 



,x 10 



t = 1 0QT 



, x 10 



t =2QQT 



,x 10 



, x 10 



1 2 

V/C „ 
s,0 

x = 1 .on - 



,x 10 



1 2 

v/c . 

s,0 

„ x=io.oa _ 

, 4 f 



I 2 

v/c „ 

s,0 

t = 40ST 



, x 10 



I 2 

v/c . 

s,0 

t =80fi~ 



,x 10 



^-ririTl I I.I I nrrnrhh-^ 




,x 10 



1 2 
V/C „ 

s,0 

x = 100.0fi~ 



2 

v/c 



s,o 




Figure 6. Histograms showing the spread of particle velocities, 
relative to the initial sound speed, for each friction time at the 
end of the t c = 10 simulation. As the friction time lengthens, the 
distribution of velocities narrows due to the diminishing effect of 
the drag force. The peak of the distribution also tends to occur 
at lower velocities. 



The peak velocity decreases with increasing cooling time, 
and the width of the distribution narrows. We therefore see 
that despite the reduction in the concentration of solids that 
arises from the increasing cooling time, the likelihood of con- 
structive collisions increases due to the decreasing velocity 
dispersion. 

As shown in figures [5] and [JJ the velocities of parti- 
cles within our simulations are typically of order the local 
sound speed in the disc. For typical disc parameters, this 
ranges from 100 — 1000ms -1 . Such high velocities are typi- 
cally associated with destructive collisions between particles 
(grain destruction is typically estimated to occur for colli- 
sions above ~ 10ms _1 ), however due to the aligning of parti- 
cle velocities within density waves the relative velocities be- 
tween particles are typically much smaller, as demonstrated 
by the small spread of velocities in figures [5] and [7] especially 
in the simulations with longer values of t c . 

This opens up the possibility of a preferred region for 
planetesimal formation at radii between 30-40AU, where the 
cooling time is low enough to allow large concentrations of 
solids to accumulate, without the collisional velocities of par- 
ticles becoming high enough to inhibit growth. 
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Figure 7. Histograms showing the spread of the Tf = 1 parti- 
cles velocities, relative to the initial sound speed, at the end of 
our simulations for the four shortest cooling times considered. As 
the cooling time lengthens, the distribution of particle velocities 
narrows, with the peak of the distribution occurring at smaller 
velocities. 
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Figure 8. Maximum particle density, relative to the mean, 
against time for a range of drift values, Av = 0.0, 0.02, 0.05 and 
0.5 for cooling time t c = 20f! -1 . Even in the presence large 
amounts of radial drift, there are still significant over-densities 
in the particle density field due to the gravitational instabilities. 



3.4 Effects of Radial Drift 

The above simulations were all conducted assuming that 
there is no radial pressure gradient supporting the disc 
(Av = in equation [8} . In reality the gas in the disc is par- 
tially supported by pressure, causing it to orbit at slightly 
sub-Keplerian velocities. The dust particles however travel 
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with approximately Keplerian velocities. The resulting ve- 
locity difference results in the particles transferring angular 
momentum to the disc and spiralling inwards towards the 
star. This radial pressure gradient could potentially negate 
the effect of the spiral density waves, causing particles to 
drift through the associated pressure maxima, reducing the 
degree to which particles are concentrated. In order to quan- 
tify the effect of this radial drift we performed several sim- 
ulations at a fixed cooling time t c , varying only the value 
of the drift parameter Aw. Figure [8] shows how the maxi- 
mum surface density of particles evolves with time for each 
value of Av considered. In each case, we see that for typical 
values of Av (of order a few percent of the sound speed), 
there is very little deviation from the case with no radial 
drift. Even for very high drift values associated with turbu- 
lent discs, the overall concentration of particles is reduced 
only slightly, suggesting that the radial pressure gradient 
does little to change the efficiency of particle concentration 
within spiral waves. 



4 DISCUSSION AND CONCLUSIONS 

In this paper we present a series of simulations modelling 
dust dynamics in self-gravitating discs using the shearing 
sheet approximation. The dust particles evolve within a 
quasi-steady gaseous spiral structure produced by a com- 
bined effect of disc self-gravity and cooling. Particles are 
coupled to the gas via a drag force with a constant frictio n 
time. As in previous related studies |Rice et al.l l2004. 2006), 
we find that spiral density waves in self-gravitating discs can 
result in significant over-densities in the solid component of 
the disc. A novel contribution of this work is, however, that 
we characterised a particle trapping capability of density 
waves for a wide range of disc cooling times and particle fric- 
tion times. Particles with friction times r/ ~ 1.0f2 -1 tend 
to be affected by spiral density waves most effectively and 
exhibit the largest and the tightest concentrations in wave 
crests. At smaller friction times, the particles closely trace 
the overall gaseous spiral structure, however exhibit lesser 
concentration in wave crests. At larger friction times, par- 
ticles are almost unaffected and decouple from the gas, so 
that they mimic the spiral structure very weakly again with 
smaller concentration in spiral waves. We also found that 
particle concentration depends on cooling time (see Fig. 4). 
The shorter cooling times (t c ^ 4057 -1 ) considered show 
similar amounts of particle concentration, as the spiral den- 
sity waves are large enough to continually trap the majority 
of the particles. For longer cooling times this is no longer the 
case, as seen in the t c — 80f2 -1 and the t c — IGOQ," 1 run, 
where the maximum particle density is lower. We expect 
this trend to continue to longer cooling times and given the 
rapid increase in t c with decreasing radius (Fig. 5.), do not 
expect disc self-gravity to play a significant role in particle 
dynamics at inner radii (r < 20AU). 

This implies that particle trapping by density waves is 
not necessarily restricte d to the very outer region s of the 
disc, a concern raised bv lClarke and Lodatol ((2009). Gener- 
ally, as mentioned above, grains with effective friction time 
Tf ~ l.Ofi -1 exhibit the largest concentration in density 
waves. When this is translated in terms of particle size we 
find that grain sizes which experience the largest density 



enhancements depend on the radius in the disc, with larger 
particles more efficiently trapped at inner radii and smaller 
at outer radii. We also analysed the dependence of particle 
velocity dispersion on friction and cooling times. Given a 
fixed cooling time, the smaller the friction time, the larger 
the average velocity of particles. Conversely, at a fixed fric- 
tion time, the larger the cooling time, the smaller the particle 
velocities are. 

Based on the trends observed in this paper, we propose 
the following potential route of grain growth and evolution 
in discs: As smaller grains can be trapped primarily at larger 
radii, this can lead to continual growth as grains secularly 
migrate inward as a result of drag forces from the gas (due to 
difference between Keplerian velocity of particles and sub- 
Keplerian velocity of gas). Small grains initially at large 
radii (~ 60AU) can become trapped in gas over-densities 
and grow through collisions with other particles. This pro- 
cess can continue until the grains reach sizes larger than the 
optimally trapped size range (t/ ~ l.Ofi ). Once particles 
reach this size, one of two outcomes may occur. If the par- 
ticles have grown sufficiently large that they are no longer 
strongly coupled to the gas, then they are no longer exposed 
to rapid inward drift, and will likely remain at that radial lo- 
cation for the remainder of the discs self-gravitating phase. 
This could result in a large number of intermediate-sized 
objects forming at mid-large radii within the disc. If how- 
ever the particles have not reached a sufficiently large size 
to entirely decouple from the gas, when particles grow to a 
point they are no longer trapped efficiently by density waves 
they will start to slowly move inward in the next optimally 
trapped size range at smaller radii, where they will undergo 
further rapid growth and so on until they reach planetesimal 
sizes. However, this process may be partially offset at large 
radii, where the effective cooling time is lower, as the rel- 
ative velocities between particles may be still large enough 
to lead to the collisional destruction of large objects. At 
inner radii (~ 20 — 30AU) the velocity dispersion is much 
lower, especially for the r/ = l.Ofi particles which are 
clumped in spiral waves most effectively. This may lead to 
a ring of material, analogous to the Kuiper belt in our own 
Solar System, building up at intermediate radii within the 
disc, where the enhancement of the solids surface density 
is significantly augmented due to the self-gravitating struc- 
ture, but the velocity dispersion is sufficiently small to avoid 
collisional destruction of growing objects. 

Another interesting question that still needs to be ad- 
dressed is the effect that large accumulations of particles 
have on the gas dynamics due to the back reaction of the par- 
ticles on gas and the effect of particle self-gravity. As we have 
shown here, particle densities in the crests of spiral waves 
can reach 100-1000 times the mean particle density. For typ- 
ical proto-planetary discs, the dust to gas mass ratio is 0.01. 
If there is sufficient mass in the size range that is strongly in- 
fluenced by structures in the gas discs, there will be regions 
where the particle density can reach values comparable to, 
or sometimes even larger than, the local gas density. Un- 
der these conditions, it is no longer appropriate to treat the 
particles as massless tracer particles. To obtain a more com- 
prehensive understanding of the particles' evolution the back 
reaction of particles on gas and the self-gravity of dust par- 
ticles must be included. For example, in the case of planetes- 
imal formation as a result of streaming instability, ultimate 
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concentration of particles in the nonlinear saturated state is 
determined by particle feedback on gas that limits explosive 
growth of particle density due to linear streaming instabil- 
ity once the local du st-to-gas ration becomes o f order unity 
in gas over-densities jjohansen fe Youdin 20q3). As regards 
the role of particle self-gravity, we saw that once the parent 
gas over-density has sheared out, particles aggregated within 
it also dissolve, whereas particle self-gravity can in principle 
keep dust clumps together and prevent them from dissolv- 
ing, facilitating further bonding of grains. For example, in 
analogous simulations of particle accumulation but in gas 
over-densities in MRI-driven turbulence, ultimate sticking, 
retainment and growth of p article aggregates occ urs in fact 
thanks to their self-gravity l| Johansen et al.ll2007h . We have 
recently run simulations that do included both the back re- 
action and the particles self-gravity, which suggest that the 
back reaction of the particles is largely negligible, and does 
not significantly change the evolution of the gas, even for 
high concentrations of dust particles where the particle den- 
sity is equal to the gas, however the particle self-gravity is 
essential in the long term evolution of the solids, with several 
large, gravitationally bound structures emerging within the 
spiral density waves. Finally, only the 2D 'shearing sheet' 
case has been considered in our analysis. Therefore we do 
not unable to observe any behaviour in the z direction, such 
as dust settling. Further work must be done to expand the 
study to the 3D 'shearing box' case to study these effects. 
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